Statistical theory of high-gain free-electron laser saturation 
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We propose a new approach, based on statistical mechanics, to predict the saturated state of a 
single-pass, high-gain free-electron laser (FEL). In analogy with the violent relaxation process in self- 
gravitating systems and in the Euler equation of 2D turbulence, the initial relaxation of the laser can 
be described by the statistical mechanics of an associated Vlasov equation. The laser field intensity 
and the electron bunching parameter reach a quasi-stationary value which is well fitted by a Vlasov 
stationary state if the number of electrons N is sufficiently large. Finite TV effects (granularity) 
finally drive the system to Boltzmann-Gibbs statistical equilibrium, but this occurs on times that 
are unphysical (i.e. excessively long undulators). All theoretical predictions are successfully tested 
by means of finite N numerical experiments. 

PACS numbers: 05.20.-y Classical statistical mechanics; 41.60.Cr Free-electron lasers; 05. 45. -a Nonlinear 
dynamics and nonlinear dynamical systems. 
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The interaction of charged particles with electromag- 
netic fields is a topic of paramount importance in a large 
variety of physical phenomena, from plasma dynamics 
to astrophysical systems. On the energy exchange be- 
tween particles and field also relies the possibility of gen- 
erating coherent and tunable radiation sources, such as 
Free- Electron Lasers (FELs). In this case a relativis- 
tic electron beam propagating through a periodic mag- 
netic field (produced by an undulator) interacts with a 
co-propagating electromagnetic wave. Lasing occurs be- 
cause the undulator field and the radiation combine to 
produce a beat wave that travels slower than the speed of 
light and can be synchronized with electrons. Among dif- 
ferent schemes, single-pass, high-gain FELs are currently 
attracting growing interests as they are promising 
sources of powerful and coherent light in the UV and X 
ranges In the high-gain regime, both the light 

intensity and the longitudinal bunching of the electron 
beam increase exponentially along the undulator, until 
they reach saturation due to nonlinear effects. Under- 
standing this saturation process is important to estimate, 
and then optimize, the performance and building costs of 
a FEL. 

Theoretical analyses usually rely on dynamical meth- 
ods in combination with detailed, but rather complicated, 
numerical simulations. In this paper, we propose a new 
approach, which is based on statistical mechanics, to 
study the saturated state of a high-gain single-pass FEL. 
We restrict our analysis to the steady-state regime, which 
amounts to neglect the variation of the electromagnetic 
wave within the electron pulse length (small electrons- 
radiation slippage). However, it is important to stress 
that because of its intrinsic flexibility, we believe that 
our statistical approach will be applicable also to alter- 



native schemes, such as harmonic generation p|. 

The starting point of our study is the Colson-Bonifacio 
model || . Under the hypotheses of one dimensional 
motion and monochromatic radiation, the steady-state 
equations for the j th electron of the beam coupled to 
radiation read: 



dz 
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where N is the number of electrons in a single radia- 
tion wavelength and z = 2k u pz^f 2 / '(7)0 is the rescaled 
longitudinal coordinate, which plays the role of time. 
Here, p = (a w Lu p /4:ck u ) 2 ^ 3 /j r is the so-called Pierce 
parameter, j r the resonant energy, (7)0 the mean en- 
ergy of the electrons at the undulator's entrance, k u the 
wave vector of the undulator, lo p — (e 2 n/me ) 1 ^ 2 the 
plasma frequency, c the speed of light, e and m respec- 
tively the charge and mass of one electron. Further, 
a w = (eB w /k u mc 2 ), where B w is the rms undulator 
field, for the case of a helical undulator. By introduc- 
ing k as the wavenumber of the FEL radiation, the phase 
is defined by = (k + k u )z — 25pk u z^ 2 / (7)0 and its 
conjugate momentum p = (7 — (7)o)/(p(7)o)- A is the 
scaled field amplitude, a complex vector, transversal to 
z, A = (A x , A y ) 7\. Finally, the detuning parameter is 
given by 5 = ((7)0 — 7, 2 )/ (2pj 2 ), and measures the aver- 
age relative deviation from the resonance condition. 

Although very simple, such a model captures the main 
features of the dynamics of the single-pass FEL, as shown 
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by a systematic comparison with numerical predictions 
based on more complete approaches 0, 0. Using this 
model, we are able to predict analytically the mean sat- 
urated laser intensity, the electron beam bunching, and 
the electrons' velocity distribution, for a wide class of 
initial conditions (i.e. energy spread, bunching and ra- 
diation intensity). The analytical results agree very well 
with numerical simulations. 

The above system of equations can be derived from the 
Hamiltonian 

N 2 N 

ff = ^|-JV« + 2VT^sin(0 j - <p), (4) 

where the intensity / and the phase (p of the wave are 
related to A = A x + iA y = \fl e~ lv . In addition to the 
"energy" H , the total momentum P = £V Pj + N AA* 
is also a conserved quantity. Let us note that one can 
always take P = 0, upon a shift in the detuning [i"o| : 
thus, we always suppose P — in the following. 

It is important to emphasize that Hamiltonian (0 
models the interaction between radiation and electrons. 
Hence, it describes a quite universal phenomenon which 
is encountered in many branches of physics. As an exam- 
ple, in the context of plasma theory, the so-called plasma- 
wave Hamiltonian (0 characterizes the self-consistent in- 
teraction between a Langmuir wave and N particles, after 
an appropriate redefinition of the variables involved |ld| . 
Establishing a formal bridge between these two areas al- 
lows to recast in the context of the single-pass LINAC 
FEL numerous results originally derived in the frame- 
work of plasma physics. In addition, Hamiltonian Q 
can be viewed as a direct generalization of mean-field 
models 0, 0, ^| , which are widely studied nowadays 
because of their intriguing features: statistical ensemble 
inequivalence, negative specific heat, dynamical stabiliza- 
tion of out-of-equilibrium structures. 

In plasma physics, it was numerically shown [Tol IhH 
that, in the region of instability, wave amplification oc- 
curs in two steps. One first observes an exponential 
growth of the wave amplitude, followed by damped oscil- 
lations around a well defined level. However, the system 
does not reach a stationary state and this initial stage is 
followed by a slow relaxation towards the final statisti- 
cal equilibrium. An example of this behavior is shown in 

Fig.m 

The separation into two distinct timescales character- 
izes also the dynamics of self-gravitating systems and 
is a well-known phenomenon in astrophysics |l5l Il6| . 
The intermediate quasi-stationary states live longer and 
longer as the number of particles N is increased. It is 
believed that galaxi es ( N ~ 10 11 ) are well described by 
Vlasov equilibrium 15], which characterizes the quasi- 
stationary state of the iV-particle system (see below). 
On the contrary, Boltzmann-Gibbs statistics applies to 
the "smaller" (N ~ 10 6 ) globular clusters. 

A typical evolution of the radiation intensity I as a 
function of the longitudinal coordinate z, according to 
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FIG. 1: Typical evolution of the radiation intensity using 
Eqs. Q-©; the detuning S is set to 0, the energy per electron 
H/N = 0.2 and N = 10 4 electrons are simulated. The inset 
presents averaged simulations on longer times for different 
values of N: 5- 10 3 (curve 1), 400 (curve 2) and 100 (curve 3). 



the Free Electron Laser model 10-10 is displayed in 
Fig. starting from a very weak radiation, the inten- 
sity grows exponentially and saturates, oscillating around 
a well defined value. This growth and first relaxation 
of the system (usually called" violent relaxation" in as- 
trophysics) is governed by the Vlasov equation 0, ^} , 
which is rigorously derived by taking the continuum limit 
(N — > oo at fixed volume and energy per particle). On 
longer timescales, whose duration strongly depends on 
the particles number N (see the inset of Fig. , there is 
a slow drift of the intensity of the beam towards the final 
asymptotic plateau determined by the Boltzmann-Gibbs 
statistics. Such process is driven by granularity, a finite- 
N effect 0, flH llrJ ] . This final relaxation takes place on 
an extremely long time scale, well beyond the physical 
constraints imposed by a reasonable undulator length. 
We thus concentrate in the following on the Vlasov de- 
scription of the dynamics. 

A linear analysis Q leads directly to the determina- 
tion of the boundaries of the instability domain, which 
are mainly controlled by the detuning S and by the initial 
energy per electron. In the case of a monoenergetic elec- 
tron beam, the instability disappears for S > S c ~ 1.9. 
Linear analysis also provides estimates of the growth rate 
of /. However, getting insights on the saturated state 
requires a nonlinear study of the system; the standard 
approach to this problem is mainly dynamical, as for in- 
stance in Ref. [lq . In the following we discuss a new 
procedure, based on statistical mechanics. 

As sketched in the previous discussion, we are inter- 
ested in the intermediate metastable state and, therefore, 
we will first consider the statistical theory of the Vlasov 
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equation,originally introduced in the astrophysical con- 
text [l5L Hq . The basic idea is to coarse-grain the micro- 
scopic one-particle distribution function f(9,p 7 t), which 
is stirred and filamcnted at smaller and smaller scales by 
the Vlasov time evolution. An entropy is then associ- 
ated to the coarse-grained distribution /, which essen- 
tially counts the number of microscopic configurations. 
Equilibrium is then computed by maximizing this en- 
tropy while imposing the dynamical constraints. A rigor- 
ous description of this procedure can be found in Ref. [ljj 
in the context of two-dimensional Euler hydrodynamics. 

In the continuum limit, Eqs. Q-© lead to the follow- 
ing Vlasov-wave system: 



df 

dz 
dA x 
dz 

dAy 

dz 



df df 
-p-^+2(A x cos9-A y sia6)-^- , (5) 



SAu 



1 

2^ 



/ cos 6 d9 dp , 



SA X - — I fsm9d0dp 



(6) 
(7) 



Note that these equations have been studied numerically 
in a recent work by Vinokurov et al. [2(j, for the case 
5 = 0. The Vlasov-wave equations (I5I7H conserve the 
pseudo-energy 



dpd6^-f(6,p)-S(Al + Al) 



+2 dpd9 f(9,p)(A x sin ( 



and the total momentum 



dpde pf(6,p)+A 2 x 



Ay COS 



0) , (8) 



(9) 



For the sake of simplicity, let us suppose that the beam 
is initially unbunched, and that energies are distributed 
according to a step function, such that 



/(0,p,t = O) = fo = jK if -P<P<P 
4irp 



= 



otherwise. 



(10) 



As far as one is dealing with small energy dispersions, the 
profile l|10(l . called waterbag initial condition, represents 
a good approximation of a more natural Gaussian initial 
distribution. Numerical tests fully confirm the validity 
of this simple observation. According to Qll)|). / takes 
only two distinct values, and coarse-graining amounts to 
perform a local average of both. The entropy per particle 
associated with the coarse-grained distribution / is then 
a mixing entropy |llL Il6| and reads 



*(/) 



dpde 



JO fo 



1 



fo 



In 1 



]_ 

fo 



As the electromagnetic radiation represents only two de- 
grees of freedom within the (2N + 2) of Hamiltonian 
its contribution to entropy can be neglected. 



The equilibrium state is computed jll| by solving the 
constrained variational problem: 

H(f,A x ,A y ) = Ne; IdOdpf = 1; 
P(f,A x ,A y )=a). (12) 



a) = max ( s(f) 

f ' ]A X ,j4y 



Introducing three Lagrange multipliers (3, A and \x for the 
energy, momentum and normalization constraints and 
differentiating Eq. ifT^ with respect to /, one gets the 
equilibrium distribution 



f = fo 



e -(i(p 2 /2+2A S \n8)-\p-,i 
I _|_ e -f3(p 2 /2+2Asine)-\p- f i, 



(13) 



By differentiating Eq. I|12|l with respect to A x and A y , 
one obtains in addition the expression for the amplitude 
of the wave 



A - x l Al + A 2 y - 



(3 



A 



dpd9 sin 9f(6,p). (14) 



Using the above equations for the three constraints, the 
statistical equilibrium calculation is now reduced to find- 
ing the values of /3, A and /i as functions of energy e 
and total momentum a. This last step, performed nu- 
merically using for example a Newton-Raphson method, 
leads directly to the estimates of the main physical pa- 
rameters. 

Furthermore, let us stress that in the limit of a vanish- 
ing energy dispersion, the area occupied by the / = fo 
level in the one-particle phase space is small, so that the 
coarse-grained distribution / verifies / <C fo everywhere. 
The second term in the entropy l|llf> is thus negligible, 
and l|13fl reduces to the Gibbs distribution 



/ oc e 



-/3(p 2 /2+2Asin0)-Ap 



(15) 



Vlasov equilibrium is in that case equivalent to the full 
statistical equilibrium. Then, solving the constraint 
equations yields 



b = A 3 ~SA 
A = I e - 6 A 2 



-A 4 
2 



6(6) , 



(16) 
(17) 



where b = \ ^2j e iej \/N is the bunching parameter and 
O is the reciprocal function of I\ (x) / Jo (x) , being the 
modified Bessel function of order n. Let us note that 
Eqs. (|16fl and i|17|) give the microcanonical solution of 
Hamiltonian The canonical solution of the same 

Hamiltonian in the context of plasma physics was ob- 
tained in [23. It turns out that the two ensembles are 
equivalent, which was not granted a priori for such a self- 
consistent system with infinite range interactions [l3j . 
Let us remark that Eqs. l(TB|) and (|T7|) were obtained 
in Ref. 0] using several hypothesis, suggested only by 
numerical simulations. Here, a statistical mechanics ap- 
proach gives a complete and self-consistent framework to 
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justify their derivation. In particular, let us emphasize 
that, contrary to the previous approach, it is not neces- 
sary to choose a priori the distribution /, which is fully 
determined by the method of solution. 
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FIG. 3: Comparison between theory (solid and long-dashed 
lines) and simulations (symbols) for a non monoenergetic 
beam when the energy e, characterizing the velocity disper- 
sion of the initial electron beam, is varied. The dotted lines 
represent the intensity and bunching predicted by the full sta- 
tistical equilibrium given by Eqs. (11611 7B . not very appropri- 
ate here, whereas the solid line and long-dashed lines refer to 
the Vlasov equilibrium defined by Eqs. 1131 and 1141 . The dis- 
crepancy between theory and numerical experiments is small 
over the whole range of explored energies. 

In the case of a non monoenergetic beam, Figure |3| 
presents the comparison between the estimates of the 
above theoretical analysis and the results of direct nu- 
merical simulations of Hamiltonian Q, after time aver- 
aging. The comparison is shown in the energy range that 
allows the amplification process to take place. The good 
agreement for intensity and bunching provides therefore 
an a posteriori, but striking, support for Vlasov statisti- 
cal equilibrium. 



FIG. 2: Comparison between theory (solid and long-dashed 
lines) and simulations (symbols) for a monoenergetic beam, 
varying the detuning 8. The vertical dotted line, <5 = 8 C ~ 1.9, 
represents the transition from the low to the high-gain regime. 

Figure |21 presents the comparison between the analyti- 
cal predictions and the numerical simulations performed 
using equations i n th e case of a monoenergetic 

beam. Numerical data are time averaged. The agree- 
ment is remarkably good for 6 < 0.5, and is accurate 
up to the threshold value <5 C , although phase space mix- 
ing is probably less effective for larger detuning. For 
8 > S c , there is no amplification, hence, both intensity 
and bunching stick to their initial vanishing values. This 
transition, purely dynamical, cannot be reproduced by 
the statistical analysis. 




In this paper, we have proposed a new approach to 
study the saturated state of the Compton Free Electron 
Laser, based on a statistical mechanics approach in the 
framework of Colson-Bonifacio's model [6(. By drawing 
analogies with the statistical theory of violent relaxation 
in astrophysics and 2D Euler turbulence, we have de- 
rived analytical estimates of the saturated intensity and 
bunching. In addition to providing a deeper insight into 
the physical behaviour of this system, the results of our 
theory agree very well with numerical simulations. Due 
to its intrinsic flexibility, it may be possible to adapt the 
statistical approach to more complete models and com- 
plex schemes, thus allowing a direct comparison between 
analytical studies and experiments on real devices. Such 
a statistical approach could be used as a tool to define 
future strategies aiming at optimizing FEL performance. 
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